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ABSTRACT 

The 2011 outburst of the black hole candidate IGR J17091-3624 followed the canonical track of state 
transitions along with the evolution of Quasi-Periodic Oscillation (QPO) frequencies before it began 
exhibiting various variability classes similar to GRS 1915+105. We use this canonical evolution of 
spectral and temporal properties to determine the mass of IGR J17091-3624, using three different 
methods, viz: Photon Index (T) - QPO frequency (v) correlation , QPO frequency (v) - Time (day) 
evolution and broadband spectral modelling based on Two Component Advective Flow. We provide a 
combined mass estimate for the source using a Naive Bayes based joint likelihood approach. This 
gives a probable mass range of 11.8 M 0 — 13.7M 0 . Considering each individual estimate and taking 
the lowermost and uppermost bounds among all the three methods, we get a mass range of 8.7 M 0 — 

15.6 Mq with 90% confidence. We discuss the possible implications of our findings in the context of 
two component accretion flow. 

Keywords: accretion, accretion disks — black hole physics — radiation mechanisms: non-thermal 
X-rays: individual (IGR J17091-3624) 


1. INTRODUCTION 

The enigmatic Galactic black hole candidate IGR 
J17091-3624 underwent multiple outbursts in 1994, 
2001, 200 3 and 2007 before the recent one in Febru¬ 
ary 2011 (lin’t Zand et all 120031 iRevnivtsev et all 120031 
iCanitanio et all 120061 120121 ). During the 2011 out¬ 
burst, the source displayed st ate trans it ions li ke ot her 
outbursting black hole so urces (jRemillard & McClintock 
120061 iNandi et al.l 120 121) from the Hard to Interme¬ 
diate states (IPahari et a.1.112011b IRodriguez et all 120111 : 
ICapitanio et all I2012t liver fe Nandil 1201317 It also 

showed presence of a state_with charact eristic s similar 
to a ‘canonical’ soft state (liver fc Nandil [20131) of other 
outbursting black hole sources. After evolving along 
the canonical track, it started to display a number of 
temporal variabili t y cla sses similar to GRS 1915+105 
(|Altamirano et ahl 120111 ) . The presence of High Fre¬ 
quency QPO (~ 66 Hz) similar to GRS 1915+105 
(jAltamirano fc Bellonill2012[l and radio dete ctions, albeit 
at the level of a few mJv (IRodriguez et all [200) in the 
Hard and the Intermediate states were also seen in this 
outburst. 

The similarity to GRS 1915+105 triggered immense 
interest about the nature of this source. As of now, 
its precise location upto 5" is known and possible op¬ 
tical _an<i infrar ed co unterparts have been identified 
dBodaghee et al.ll2012f ) within the Chandra error circle. 
The distance to the source, mass and orbital parame¬ 
ters of both objects in the binary are currently unknown. 
What is certainly known, is that IGR J17091-3624is only 
the second black hole candidate to display a wide range 
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of temporal and spectral variations after GRS 1915+105, 
althou gh at much lower (10 t o 50 times) observed flux 
levels jAltamirano et aD 1201 111 . More knowledge about 
IGR ,117091-3624 should allow us then to make some 
meaningful comparisons between these two sources. Mul¬ 
tiple attempts have been made before to determine the 
mass of IGR J17091-3624. T hese attempts have pla ced 
it anywhe re between < 3 M^ (jAltamirano et al. 201111 to 
< 15M 0 (lAltamirano fe Bellonill2012l) . Other attempts 
to es timate the mass are < 5 M^ (iRao fc Vadawalel 
l2fiT2l). 9.13 ± 2.25 M .-. (IPahari et al.l 1201 111 , and ~6M^ 
(iRebusco et al.ll2012l) . The wide range of these estimates 
doesn’t enable us to even pin down whether the source 
is a massive stellar black hole binary or an exceptionally 
low mass black hole candidate. 

In this paper, we attempt to make an estimate 
of the mass with tighter constraints using X-ray ob¬ 
servations made during the rising phase of its 2011 
outburst. Here, we discuss three different methods 
(two of which are independent of distance to the 
source) to estimate its mass. In the 3 rd approach, we 
study the broadband (0.5 to 100.0 keV) spectrum us- 
ing t h e two com ponent a dvect i ve flow (TCAF) m odel 
(IChakrabarti fc Titarchuklll995l : iChakrabarti fe Mandall 
I2006H . We use these three estimates to put a limit on 
the mass of the central object. Finally, we also discuss 
the results obtained from our spectro-temporal analy¬ 
sis in the light of two different types of a ccretio n flows 
i.e., a Keplerian disk (iShakura fe Sunvaevlll973l l on the 
equatorial pl ane, sand wiched b y a su b-Kepleria n flow 

jlChakrabar^&^itarchu^U^iJSIl&kskSiii^^^il^J 

20061: IWu et al. I 120021 ISmith et al.l I2001L 120021 l2007t 
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disk in the vicinity of an accreting black hole. 

2. OBSERVATIONS AND ANALYSIS 

We use archival data of observations made by the 
RXTE , Swift and INTEGRAL satellites for the 2011 out¬ 
burst (from 06 Feb (MJD 55598.28) onwards). We re¬ 
strict ourselves to the rising phase of the outburst, before 
the source entered into its enigmatic and unpredictable 
variability phase. In all, we analyze ~ 40 days of data 
from these observatories. 

2.1. Swift data reduction 

We use the XRT data from the Swift suite of instru¬ 
ments. The XRT data consists of windowed timing 
mode observations, which reduces pile-up effects from 
the source. The maximum observed count rate is about 
40 counts s" 1 2 3 4 . This ensures th at the XRT data is not 
piled up ([Romano et, al. I [2 00 fill . We analyse these XRT 
data-sets by following the steps as given in the XRT man- 
uafl For this purpose, we use HEASoft 6.15.1 and its 
associated ftools packages. We restrict ourselves to use 
events of grade 0-2, while extracting timing and spec¬ 
tral data. The source and background regions are ex¬ 
tracted by selecting a circular region of diameter 40 pix¬ 
els (~ 90") from the data images. While doing spectral 
modelling, we found that the low energy spectrum (0.5 
-1.0 keV) of the Swift XRT data did not fit well and 
showed some excess in the soft state. This was due to 
uncertainty in position of the source in the binned pixels 
of WT mod(0. By using position dependent rm/fl made 
for different probable positions of the source on the de¬ 
tector Y pixel, we accounted for this apparent excess in 
the soft spectrum (K. Page 2014, private communica¬ 
tion). We use such position dependent rmfs for each of 
our Swift XRT data-sets. 

2.2. INTEGRAL data reduction 

We use the IBIS/ISGRI instrument dLebrun et al. I 
I2003T) from the INTEGRAL suite. We extract spectral 
information from ISGRI by following the steps as given 
in their data reduction manuaQ. To obtain the spec¬ 
trum, we use OS A v 10.0 along with updated versions 
of the calibration and response files. The spectrum from 
ISGRI is rebinned in order to get more data points for 
modelling the spectrum (see Table [T|). We extract spec¬ 
tral information from ISGRI only for those dates which 
have simultaneous Swift observations. This helps us to 
obtain broadband (0.5 - 100 keV) data for our spectral 
analysis. In all the data-sets that we use, IGR J17091- 
3624 is visible in the ISGRI detector with significance 
greater than 7a. However, owing to the low source flux, 
there is non-detection or very low significance detection 
in the JEM-X detector. Hence, we do not use any spec¬ 
tral data from JEM-X. 

2.3. RXTE data reduction 

We use the PCU2 and HEXTE d etector data of the 
RXTE satellite. However, as noted (lPahari et all 120111 

1 Swift XRT Data Reduction Guide, vl.2, 2005 

2 http://www.swift.ac.uk/analysis/xrt/digest_cal.php+abs 

3 http://www.swift.ac.uk/analysis/xrt/rmfs.php 

4 IBIS Analysis User Manual, Chernyakova, M. et al., 2012, 
ISDC 


iCapitanio et all l2012t liver fe Nandil 120131) , the observa¬ 
tions made by RXTE before 23 Feb are contaminated 
due to the nearby LMXB NS source GX 349+2. We 
do use these contaminated observations for our timing 
analysis. Specifically, we use them to find the frequency 
of QPOs from the power density spectrum (PDS). The 
source GX 349+2 does not exhibit any QPO like fea¬ 
tures in its PDS. We investigated the observations car¬ 
ried out by the JEM-X payload onboard INTEGRAL 
on 11 Feb 2011 to ascertain this. Also, as no t ed by 
IQ’Neill et all (120021 1: lAgrawal fe Bhattaclrarvval ([20031 1 
no QPOs have ever been found in GX 349+2, in-spite of 
extensive searches carried out over all the spectral and 
temporal states of this LMXB. This leads us to believe 
that the QPO like features in the PDS obtained from 
PCU2 data are indeed from IGR J17091-3624. We use 
only the non contaminated data from PCU2 and HEXTE 
detector for our spectral analysis. 

2.4. Temporal Analysis 

The Power Density Spectrum (PDS) of all data-sets is 
computed in the energy bands afforded by the observ¬ 
ing instruments, viz. Swift XRT (0.5 - 10.0 keV) and 
RXTE PCA (3.0 - 30.0 keV). The PDS are obtained us¬ 
ing GHATS vl. 10. a customized IDL based timing pack¬ 
age, which takes care of re-binning, Poisson n oise su b- 
traction and dead-time correction ([Zhang et al. lj 19951 ) of 
the FFT data. Each of the PDS is made by using tem¬ 
poral data of bin size 3.52 ms (for Swift XRT data) and 
0.244 ms (for RXTE PCA data). All observations are 
divided into segments of 128.0 s and a PDS is computed 
individually for each such segment. We then obtain the 
final PDS for an entire observation by averaging the PDS 
from each individual segment. Finally, we systematically 
search each averaged PDS for the presence of QPOs. A 
combination of Lorentzian features is used to model the 
PDS. We select only those QPO features which have a 
significance greater than 3er for our analysis. The Q- 
factor (v/FWHM) of these selected QPOs is found to 
vary between 3 to 10. The rms power under the QPO 
varies in between 8% to 17% and the QPO frequencies are 
observed to increase with time from 0.055 Hz to ~ 5 Hz 
(lSha,oo.shriiko\l[20TTl lPahari et alJfeOllllRodriguez et all 

12011b liver fc Nandil l2013[ h We have modelled this evo¬ 
lution of QPO frequencies in 33.21 The error on QPO 
frequency is taken to be ler deviation from its centroid 
frequency, which is obtained by scaling its FWHM. 

2.5. Spectral Analysis 

The spectral analysis is carried out using XSPEC 
vl2.8. lg. We do spectral modelling for individual obser¬ 
vations of each instrument. Apart from this we also anal¬ 
yse broadband spectra (0.5 keV to 100.0 keV) as obtained 
by choosing simultaneous observations from Swift XRT 
and INTEGRAL IBIS. We rebin the spectrum obtained 
from each instrument and add systematic errors as men¬ 
tioned in Table Q] before fitting the data-sets. We model 
the energy s pectrum by two different methods. The first 
method (see ICapitanio et al.l 120121 liver fc Nandil 1201311 
is a phenomenological model fitting of the spectrum and 
consists of a diskbb and a powerlaw (or cutof fpl) com¬ 
ponent. In this phenomenological modelling, we find that 

5 http: //www.brera.inaf.it/utenti/belloni/GHATS_Package/Home.html 


































Determination of mass of IGR J17091-3624 


3 


the powerlaw photon index (r) increases with time and 
is positively correlated with the QPO frequencies ob¬ 
tained in 1 12.41 We model this correlation as presented 
in 93.11 to obtain the source mass. We use ler errors 
on fit parameters as estimated from the XSPEC error 
command for our model fitting routines. In the second 
method of spectral modelling, we calculate a model spec¬ 
trum for each broadband observation using Two Com¬ 
ponent Advective Flow (TCAF). This is done by in¬ 
cluding radiative hydrodynamics of the accreting plasma 
self-consistently into the go vern ing equ ations of the flo w 
(IChakrabarti fc TitarchuHll995l : IChakrabarti fc Mandall 
l2006h . We then use this calculated spect rum to fit the 
observed spectrum, which is presented in 93.31 

3. MODELS AND RESULTS 

We consider three different approaches to estimate 
mass of the central source. These methods are discussed 
in the following sub-sections: spectro-temporal corre¬ 
lation 1 93.ip . the evolution of QP O fre quencies 1 93.21) 
and broadband spectral modelling ( 93.31) . Then, we use 
these three different estimates to derive a single set of 
bounds for the m ass. The estimates a re m a de u nder the 
TCAF paradigm (IChakrabarti fe TitarchuUl995f) , which 
enables us to model both evolution of spectral and tem¬ 
poral features under a single paradigm, thereby ensuring 
consistency of the estimate. In this paper, we have used 
radial distance in units of Schwarzscliild radius (r g ) and 
mass in units of solar mass (M 0 ). 

3.1. Photon Index (T) - QPO frequency (v) correlation 

A model for the observed correlation between the spec¬ 
tral fitted photon index (r) and the observed QPO fre¬ 
quency (v) for black hole cand i dates is p resented in 
ITitarchuk fc Fioritol (120041 1. iShaposhnikov fe TitarchuS 
(l2007t ) (hereafter ST07) have presented an implicit em¬ 
pirical scaling relation between the correlation curves of 
different black hole sources and used it to estimate the 
mass of a few such sources. The model uses a ‘Compton 
cloud’ around the central source and a transition layer 
/ region between this Compton cloud and the Keple- 
rian disk. Any change in the size of the cloud leads to 
a change in both the QPO frequency (explained as the 
magneto-acousti c reso na nce osc illations of the bounded 
transition layer (ITitarchuk fe Woodll2002lf l and the pho¬ 
ton index of the emitted spectrum (T) (due to vary¬ 
ing optical depth of the hot electrons). We believe 
that the ‘Compton cloud’ in this model is the same 
as the post-shock region of t he accretio n disk unde r 
the TCAF paradigm (|Cliakrabarti fe Titarchukl 119951) . 
The TCAF paradigm too demonstrates both the steep¬ 
ening of spectral i ndex and the increase in QPO fre¬ 
quen ci es with time (Chakrabarti et al. 2008; Nan di et al.1 
l2012t iRadhika fo Nandil 12014 iDebnath et al. I I2014H . 
In addition, t he accret ion disk co nfig u ratio n used in 
ITitarchuk fc Woo~dl (120021) ; ITitarchuk fc Fioritol (I2004T ) is 
similar to that under the TCAF paradigm. We thus ex¬ 
pect the empirical correlation of ST07 should hold under 
the TCAF paradigm as well. 

The correlation is fitted using the empirical relation 
given in ST07 : 


T - QPO Correlation 



Figure 1. Photon index (r) - QPO frequency (z/) correlation plot. 
Diamond points (blue) are Swift XRT QPOs. Rest of the QPO fre¬ 
quencies are from RXTE PC A. Bottom panel shows the weighted 
residuals in each point. The weights are assigned as given in Ap¬ 
pendix I. The best fit curve (red solid line) using Eqn. [T]gives mass 
to be 10.90“( g® Mg. 

where A is the value at which the photon index sat¬ 
urates, v tr is the threshold frequency above which the 
levelling off/saturation of T happens and B is the slope 
of the correlation which scales with mass. The parame¬ 
ter D controls how fast (i.e., over what frequency range) 
the transition occurs. We note that this is an empiri¬ 
cal scaling equation relating the shape of the correlation 
curve with the masses of different black hole systems. 
To find the mass of an unknown source requires com¬ 
parison of its correlation curve with a reference curve of 
a source with known mass. An inherent assumption in 
this scaling equation is the similarity of the shape of the 
corre lation curves between t hes e two sources (see also 
IShaposhnikov fc Titarchukll200 9 ! f. The slope of this cor¬ 
relation curve, as the source transits from low photon 
index and frequency values in the harder states to the 
higher saturated photon index and frequency values of 
the softer states, is related to the mass of the source. 
The empirical scaling relation estimates this slope and 
a comparison of slopes from different sources gives a di¬ 
rect comparison of their masses. We consider the ris¬ 
ing phase of the 2005 outburst of GRO J1655-40 as the 
reference correlation curve for estimating mass of IGR 
J17091-3624. This correlation was successfully used in 
ST07 to predict the mass of GRS 1915+105, a source 
with a very similar correlation curve to IGR J17091-3624. 
We keep the value of D fixed at 1.0 as done in ST07. We 
note though, that changing the value of D changes the 
final value of mass by a small amount. The best fit¬ 
ted curve (for D = 1.0) using Eqn. [l] is shown in Fig. [T] 
along-with the residuals (bottom panel). This fitting has 
been do ne using C raig Markwardt’s IDL based routines^ 
(lMarkwardtll2009D , suitably modified for accommodating- 
errors in both variables, as is explained further in Ap¬ 
pendix I. This fitting gives us values as A = 2.38^'°^, B 
= 0.23^qq3 and v tr = 4.47^gg7 Hz. The best fit mass (as 
obtained by scaling the B parameter with GRO J1655- 
40’s mass which is taken to be 6.3 ± 0.5 M 0 (I Green et al.l 
120011) ) comes as Mbh = 1O-9 O+u 67^0 with a fit statistic 


r(i/) = A- DB In 


exp \ Vtr D 1 | + 1 


(1) 


http://purl.com/net/mpfit 
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Table 1 

Reconditioning of reduced spectral data-sets 


Instrument 

Binning 

Systematic error 

Swift XRT 

Minimum 20 counts per bin 

3% 

INTEGRAL IBIS/ISGRI 

25 bins between 13 keV - 150 keV“ 

5% 

RXTE PGA 

None 

0.5% 


a IBIS/ISGRI has Gaussian distributed errors and doesn’t need rebinning for 
noise statistics. By default 11 bins are present between 13 keV - 150 keV. 


(Appendix I) of O = 16.964 for 18 degrees of freedom. 

3.2. QPO frequency (v) - Time (day) evolution 

We model the monotonic increase in QPO frequen¬ 
cies with time using the Pro pagatin g Oscillatory Shock 
(POS) model (IC hakrabarti &; Manickamll2000t hereafter 
CM00; IChakrabarti et all 120081 ). This model which 
comes under the TCAF paradigm, explains formation 
of QPOs due to an oscillating shock front formed in 
the sub-Keplerian component of the flow. The shock 
front starts to oscillate at infall time scales, where the 
infall time is the time taken by material parcel to fall 
onto the central source from the shock location. The 
reasons for shock oscillation can be either resonance 
oscillations, when the infall time scale is comparable 
to the coolin g time sc ale of the Comptonizing post¬ 
shock region dMolteni et al.lH996h . or unstable pertur- 
bations occurring in the shocked flow due to its viscosity 
(|Lee et al.ll2011l : lDas et al.l 120141) . The frequency of the 
QPO, which is the rate at which the shock front oscil¬ 
lates (fChakraba^^^^^nickam 200( 1 ICh akrabar ti et all 
120081120091 : iNandi et al.ll2012URadhika fc Nandill2014H is 

given as : 



55596.0 55603.2 55610.5 55617.8 55625.0 

MJD (days) 


Days (tO = 06 Feb 2011) 

-3 2 8 14 20 26 



55597.0 55604.0 55611.0 55618.0 55625.0 


_ VsO _ 

27 jRr s y/r s — 1.0 


( 2 ) 


where, r s is the location of the drifting shock front, v s $ 
is the inverse of light-travel time across the black hole, 
given as v s $ = — and R is the compression ratio (i.e., 

ratio of post-shock to pre-shock densities of the flow). 
In Eqn. [2j we have introduced an additional 27 t factor 
in the denominator to ensure that the entire axisymmet- 
ric shock front region (27rr s ) is oscillating in-phase and 
thi s is consist ent with the hydrodynamic simulations of 
iMolteni et al.l (Il996i l. The shock front drifts inward with 
a constant velocity vq giving, 

r 3 (t) = r s0 ± v 0 (-——); r s > r s , mm , (3) 

r 9 

r s (t) =r s ,min] otherwise, (4) 


where r s o is the initial shock location at to (1 st day 
of the outburst). We think that the observed flatten¬ 
ing/saturation of QPO frequencies (see Figure [2j may 
be due to the fact that the oscillating shock front stops 
drifting further inward. We account for this in the equa¬ 
tion by using a minimum shock location r s ^ m i n . 

We do the QPO evolution fitting by again using 
Markwardt’s IDL routined^. The difference between 
our attempt and previous studies of the POS model 
(|Chakrabarti et al.ll2008l 120091 : INandi et al.ll2012D is that 


MJD (days) 

Figure 2. Evolution of QPO frequencies with time (days) of the 
onset-phase of the 2011 outburst of IGR J17091-3624 is fitted with 
POS model solution (top panel). Marked points (with cross) in top 
panel are Swift XRT QPOs. The corresponding size of the oscillat¬ 
ing region is shown in the bottom panel with diamonds (blue) show¬ 
ing t he ‘C ompton cloud’ size as determined from spectral modelling 
(see COli . The best fit curve (red solid line) of Eqn. 2 to 4 gives 
an estimate of mass to be 14.37(j(g y( Mq. 

we have used mass as a free parameter. Most of the pub¬ 
lished work before this has dealt with establishing the 
POS model for sources of known mass by demonstrating 
that the QPO evolution can indeed be explained using 
the simple parametric expression in Eqn. [2] In our cur¬ 
rent implementation, the Schwarzschild radius r g which 
is proportional to the mass of the source is used in the 
fitting function (Eqn. OH), as a free parameter. Hence 
mass is calculated by determining r g . The other free pa¬ 
rameters in the fitting are r s q (the initial shock location) 
and rs^min- The first day of the outburst (to) is set from 
observations at 06 Feb 2011 (MJD 55598.28). R is held 
fixed at 3.0 and vo at 10 m/s, which are typical values 
for these parameters. This is inferred from the fact that 
the shock front oscillat es to give a detect able QPO for 
R between 2.0 and 4.0 (INandi et al.ll2012h . The typical 
inward drift velocity of this front as obta i ned from pre¬ 
vious studies lies between 5 - 15 m/s (ICha krabart i et all 
12001 l2009t INandi et, al.l I Ml IRadhika fc Nandill20l1h 
We use the mean of these ranges in our fitting and inves- 


















































Determination of mass of IGR J17091-3624 


5 


tigate the effects of variation of R and vq in the discussion 
section. 

The results of the POS model fitting is shown in Fig¬ 
ure [2] along with residuals (bottom panel of the top 
figure). The bottom figure shows the variation of the 
Comptonizing region (i.e., the size of the oscillating re¬ 
gion) which decreases as the QPO frequency increases 
with time. The dashed dotted vertical line marks the 
t ransit ion of hard-intermediate to soft-intermediate state 
(liver fc Nandil 12013 ). where photon index saturates over 
QPO frequency (see Fig. [2] and Fig. [[]). We obtain the 
fitted value of the mass to be Mbh = 14.37^g;®[ M 0 with 
Xred = T03 ( X 2 /dof = 15.43/15). 

3.3. Spectral modelling based on TCAF 

The spectral modelling is based on steady state hydro- 
dynamic solutions of the equations governing the flow of 
accreting material under the TCAF paradigm. The so¬ 
lution for a set of free parameters gives a model spec¬ 
trum, which we use to compare and fit the observed 
spect rum. We have incorpor ated the TCAF m odel 
(ICha krabarti fe Titarch ukl 119951: [Chakrabarti fc Mandall 
120061 : iM andal fc Chakrabartill2010l ) in XSPEC as a local 
model~ lArnaudlll996i l to achieve this spectral fitting. In 
this model, we consider a sub-Keplerian flow of accretion 
rate rrih on the top of a Keplerian flow of accretion rate 
rhd at the equatorial plane. Both the accretion rates 
are measured in units of the Eddington accretion rate 
(triEdd)- The other free parameters are mass of the cen¬ 
tral object (Mbh), shock location (r s ) and compression 
ration (R). A similar implementation in validating the 
TCAF mod el for a few blac k hole sources has been car¬ 
ried out bv iDebnath et al.~ (12014} ). But these previous 
studies are limited within 3.0 keV - 30.0 keV whereas 
in the current work we model the broadband spectrum 
between 0.5 keV -100 keV. 

To implement the local model in XSPEC, we use the 
HEASoft facility to create a table model by generating 
the model spectra for a large range of input parameters 
and saving it in an array. This table model can be run 
from XSPEC (XSPEC User’s Guide) to calculate the fitting 
parameters. However, the values of the fit parameters 
so obtained may not be very accurate as XSPEC does an 
interpolation to fit the data. Hence, we have done the 
spectral modelling in two steps. First, we use the table 
model to constrain the parameters in reasonably nar¬ 
row interval and then we run the actual hydrodynamic 
source code from XSPEC with the parameters previously 
obtained narrow interval to get the best fit. 

During the 2011 outburst, the source evolved from 
hard to soft state via intermediate states. We have fit¬ 
ted the broadband spectra (0.5 - 100.0 keV) for all the 
observations where we could find simultaneous or quasi- 
simultaneous broadband observations. The observations 
so found lie in each of these spectral states: hard state 
(HS) observed on 08 Feb, hard intermediate (HIMS) on 
20 Feb, Intermediate state near transition from HIMS to 
SIMS on 22 Feb, soft intermediate (SIMS) on 28 Feb, 
Intermediate state near transition from SIMS to SS on 
08 M ar and soft state (SS) on 12 Mar (see liver fe Nandil 
12011 . We model these six spectra by using our local 
TCAF model alongwtih phabs model of XSPEC to account 
for the interstellar absorption. In the model spectrum 


computation, the soft photons are supplied by the Kep¬ 
lerian disk and are inverse-Comptonized by hot electrons 
in the post-shock region. Hence, the resultant spectrum 
is a combination of soft and the hard components and 
the relative normalization is determined by the fraction 
of photons intercepted by the post-shock region. We fi¬ 
nally fit the resultant spectra with observed spectra as 
shown in Fig. 3. The TCAF model fitted parameters 
are presented in Table [2] We have kept the compres¬ 
sion ratio fixed (R=3) for all spectral fitting. We see 
that as the source moves to soft state, the shock front 
moves inward and increases while rhh decreases ex¬ 
cept for observations near the state transition bound¬ 
aries (see )2] for details). This means that in the soft 
state, a lot more soft photons come from Keplerian disk 
and there exist lesser number of hot electrons to cool 
down. This is consistent with the physical picture of 


an at 

1995 

:creting black hole source (IChakrabarti & Titarchuk 
: Smith et al.l 20011 2002), Chakrabarti & Mandal 

2006 

: Mandal & Chakrabarti! 20101). 


In our spectral fitting, we find that the column den¬ 
sity (nn) varies between (0.97 — 1.7) x 10 22 cm -2 while 
the source transits from the hard to the soft state. A 
higher value of nn in t he soft state m ay be due to en¬ 
hanced disk winds (see lKing et al .1120121) or other sources 
of opacity. In our present work, we do not assume any 
distance to the source. We calculate the overall normal¬ 
ization parameter (fixed for all spectral fits) from our 
spectral modelling. For every data-set, we calculate the 
normalization constant associated with the best fit and 
then take an average over all normalizations from indi¬ 
vidual data-sets. Finally, we fit all the data-sets again 
with this new averaged normalization constant. Also, 
we do check that the final estimation of mass from in¬ 
dividual data-sets with fixed overall normalization lies 
within the range of mass as determined by the individ¬ 
ual normalization before taking the average. The overall 
normalization depends on the distance to the source and 
inclination angle of the system. The reported distance 
to the so urce and its inclination angle varies between 10 


- 20 kpc (Altamirano et al. 

2011tlRao & Vadawale 2013T) 

and 50° - 70° (King et all 

2012) respectively. Hence for 


an inclination angle in the range (50° - 70°), the final 
normalisation constant translates into a distance range 
(10.6 - 14.6) kpc. For lower values of inclination angle, 
the distance to the source reduces even further. This sug¬ 
gests that the source flux is low, not due to the source 
being located at a very large distance. The spectral 
fitting gives a mass estimate as listed in the last col¬ 
umn of Tabled with x 2 ed (y 2 /dof) for different states 
as 160.57/157 (HS), 132.85/158 (HIMS), 238.98/269 
(HIMS-SIMS), 117.53/138 (SIMS), 341.85/283 (SIMS- 
SS) and 163.83/176 (SS) respectively. 

4. DISCUSSION AND CONCLUSIONS 

In this paper, we have tried to constrain the mass 
of IGR J17091-3624 by attempting to explain the X- 
ray spectral and timing observations under the TCAF 
paradigm. We find that the source behaviour can be 
consistently explained if the source mass is > 10 M 0 . To 
quantify this further, we first examine the limitations on 
each of our methods one by one. Finally, we attempt to 
put a single set of limits on the source mass from the dif¬ 
ferent estimates presented in this paper. The summary 
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Figure 3. Broadband energy spectra fitted by TCAF model for different states mentioned on the respective figures (residuals are shown 
in bottom panels). HS, HIMS and SIMS spectra are taken by using simultaneous data from Swift XRT and INTEGRAL IBIS. SS spectrum 
is taken by using simultaneous data from Swift XRT and RXTE PCA. The fit parameters along with the estimated mass are presented in 
Table [2] 


Table 2 

Results of the fitting for Mass estimation 


Method 


Details and Assumptions 

Mass estimate 

(A) T-u Correlation 


See ST07 

10.90+i;« Mq 

(B) QPO Evolution (POS model) 

See CM00, Nandi et al. (2012) 

14.37 + 0.67 Mq 

(C) Spectral Modelling (TCAF) 

State 

Shock location 

Accretion rate 

Mass 



(r s ) 

[m h , m d ) 

(Mq) 


HS 

314 ± 44 

0.66 + 0.04, 0.08 + 0.01 

12.7 + 0.6 


HIMS 

77.3 + 9.3 

0.58 + 0.04, 0.12 + 0.02 

10.7 + 0.5 


HIMS-SIMS 

35.1 ± 11.3 

0.16 + 0.015, 0.21 + 0.14 

11.5 ± 1.5 


SIMS 

28.2 + 4.5 

0.28 + 0.05, 0.63 + 0.08 

12.9 + 1.1 


SIMS-SS 

24.0 + 11.6 

0.12 + 0.02, 1.09 + 0.09 

13.1 ± 1.1 


SS 

19.0 + 3.7 

0.25 + 0.06, 0.64 + 0.12 

13.1 ± 1.6 

Final Bounded Mass 

Joint Likelihood : 11.8 Mq - 13.7 Mq 


of mass estimates by Photon Index (T) - QPO frequency 
(v) correlation, QPO frequency (v) - Time (day) evolu¬ 
tion and Spectral modelling based on TCAF is given in 
Table [2] 

In the T-QPO correlation model, we have noted that 
our inferred value of mass depends on the value of ‘D’ 
parameter. In this method, we find that the range of 
frequency over which the transition (bottom left to top 
right) in the correlation curve (see Fig. [lj occurs for the 
reference source of our choice, GRO J1655-40 (~ 0.1 Hz 
to 15 Hz), is about three times the range for IGR J17091- 


3624 (~ 0.05 Hz to 5 Hz). This would give us an estimate 
of D=0.33, obtained as a ratio of the frequency range of 
the two sources, as D controls the range of frequencies 
over which the correlation curve shows the above men¬ 
tioned transition. However, we also note that the value 
of D = 1.0 is used to scale the mass of GRO J16555-40 
to GRS 1915+105 (see ST07 for details). GRS 1915+105 
has a similar span of QPO frequencies like IGR J17091- 
3624. To resolve this, we redo the fit for different values 
of parameter D from 0.33 to 1.0, and note the variation 
in mass as the uncertainty / limitation of this model. 
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We find that the mass estimate shifts from a minimum 
at 9.50 Mq to a maximum of 10.90 Mg. The value of 
1.4 Mg (10.9 Mq - 9.5 Mg) is total uncertainty on the 
estimated mass using this method. To use this uncer¬ 
tainty in our final estimate, we take the lcr uncertainty 
to be one third of this value (i.e., we take 1.4 Mg to 
encompass 3<r (or 99%) of the total uncertainty of this 
method). The scaling law used to o btain t he mass from 
this method is based on the model of lTitarchuk fe Fioritol 
( 2004 1 and the empirical relation of ST07. If this empiri¬ 
cal relation and the subsequent linear scaling mis-specify 
the actual physical relation involved, then that could in- 
trod uce anoth er s ource of uncertainty in this estimate 
(see lYang et al. I (1201211 to get an idea of model uncer¬ 
tainties due to mis-specified scaling assumptions). We 
do not explicitly account for this in our paper, as further 
investigation of this model would require cross validat¬ 
ing the model against multipl e black hole system s with 
a large sample of data fsee IKessler et al.l (12013H for an 
example of uncertainty estimation). 

In the QPO evolution model, we see that the the range 
of values over which R and vo can potentially vary adds 
a source of uncertainty to our quoted values of mass. As 
done above, we try to estimate this by redoing the QPO 
evolution fit for different R (from 2.0 to 4.0) and vq (from 
5 m/s to 15 m/s). We find that this gives an uncertainty 
in mass of 1.7 Mg from the mean value (i.e., the standard 
deviation (lcr) in the final set of mass values is 1.7 Mg). 
Secondly, and more importantly, the data does not fit 
well if vq is less than 8 m /s or great er than 12 m/s. 
In a recent paper, iMondol et al.l (120151) have mentioned 
that the velocity of propagation of the shock front can 
be calculated from spectral modelling by estimating the 
shock location. However, in the present paper we treat 
spectral modelling and QPO frequency modelling as two 
independent methods and hence do not use the results 
from one method in the other method. 

The spectral fitting method under TCAF paradigm 
also has a few limitations as listed below. Our 
spectral modelling includes black body and inverse- 
Comptonization components. These two components are 
good enough to fit the spectra of IGR J17091-3624. But 
there are other sources where these two components are 
not sufficient and we need to include components like 
a Gaussian emission line profile (e.g. GRS 1915+105), 
and absorption smear edges (e.g. XTE J1859+226). We 
need to include other physical processes in our model to 
generate such components. The hydrodynamic solutions 
used to calculate the model spectra are not self-consistent 
transonic solutions. The number of fit parameters can be 
reduced further if we use transonic solutions. To account 
for the shock transition from pre-shock to post shock re¬ 
gion, we fix the value of R = 3.0 (as done in T3.2I) . We 
find that spectral modelling is not very sensitive to R in 
our case and do not account for it separately. Secondly, 
an uncertainty in overall normalization (fixed for all data¬ 
sets) may lead to the uncertainty in the estimated mass. 
Our spectral modelling constrains the system in such a 
manner, that we address the spectral shape along with 
the overall luminosity simultaneously for multiple data¬ 
sets. In principle, we do not expect the normalization 
to change across data-sets. Hence, we expect this uncer¬ 
tainty to be small enough to not significantly affect the 


mass estimation. 



Figure 4. Combining the estimates. Figure shows the results 
from the joint likelihood based approach for combining the mass 
estimates. See Appendix II and Appendix III for details. 

Given these three estimates and considering the lim¬ 
itations of the methods, we try to combine the results 
to get the overall range of mass of the object. For 
combining the estimates, we calculate the probability 
distribution function (PDF) of the mass of the central 
object from the least-squares fit for each method (see 
Appendix II). We combine results from the indepen¬ 
dent methods only. Thus, we do not combine the re¬ 
sults from r — v correlation and from QPO-day varia¬ 
tion, since both these methods share the common vari¬ 
able v (the QPO frequency). Instead we combine re¬ 
sults from §3.1 with results from §3.3, and results from 
§3.2 with results from §3.3 separately. Finally, we quote 
the lowermost and uppermost values from such a com¬ 
bination as the final mass bound. The approach that 
we take to combine the estimates is based on estimation 
of joint likelihood using the Naive Bayes theorem from 
independent PDFs (see Appendix III for details). We 
present the combined estimate in Figure U where red, 
blue and green lines are individual PDFs from §3.1, §3.2 
and §3.3 respectively. The combined estimates are repre¬ 
sented by shaded regions. The combined estimate using 
this approach for §3.1— §3.3 gives us a mass range of 
11.8 M 0 — 13.1 Mg and for §3.2— §3.3 gives us a mass 
range of 12.5 Mg — 13.7 Mg. This gives us a limit on the 
black hole mass of 11.8 Mg — 13.7 Mg. For a worst case 
estimate, we find the 90% confidence value limits for each 
estimate from its respective PDF. The lowermost and up¬ 
permost of these limits places the mass between 8.7 Mg 
- 15.6 Mg. 

Whichever way we combine the estimates, we see that 
the mass is greater than 8.7 Mg. This suggests that 
IGR J17091 also harbours a black hole with mass s imi¬ 
lar to the black hole of GRS 1915+105 fsee iReid et al.l 
120141 . for a recent mass estimate). The reason for any 
distinction between these two sources, may then be the 
difference in mass accretion rates. For IGR J17091-3624, 
we obtain the total mass accretion rate to be in the range 
(0.4 - 1.1) times the Eddington accretion rate (see Ta¬ 
ble [2]). We obtain near-Eddington accretion rate close 
to the transition from SIMS to Soft State and lowest 
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accretion rate near the transition from HIMS to SIMS 
(liver fc Nandill2013li . We note that the non-simultaneous 
nature of the spectra in these two observations with a gap 
of ~ 6 hours may have contributed to this. This may also 
be due to many other time depen dent n onli near effect s 
(for e.g. flow viscosity (iManda fc Chakrabartil [2010b 
or jet activity (Radhik a fc Nandi 120141) ') involved in the 
system during the state transitions. Our spectral fit¬ 
ting is done by using a steady state model and does not 
account for these time varying phenomena. We have at¬ 
tempted to apply our spectral modelling techniques on 
GRS 1915+105 as well for a relative comparison of flow 
accretion rate. While modelling the x class spectrum of 
GRS 1915+105, we find evidence for super-Eddington 
accretion rates (~ 6.5 times the Eddington accretion 
rate) in the system. In future, we would like to do a 
detailed broadband spectral modelling of different vari¬ 
ability classes of GRS 1915+105 and IGR J17091-3624to 
make a comparative study of the variation in accretion 
rate in both systems. However, it is suggestive to note 
that the difference in accretion rates, could be the reason 
for the observed difference in X-ray flux. 

If that indeed is the case, then it raises the question 
of how both sub- and super-Eddington accretion rates 
give rise to similar kinds of multiple variability classes 
as seen in these objects. It is possible that the com¬ 
plex variability classes observed in these sources (IGR 
J17091-3624 & GRS 1915+105) could be due to the in¬ 
terplay between the Keplerian and sub-Keplerian flows in 
p resence of ou tflo ws/ winds in the soft-inter mediate state 
IChakra barti fc Titarch ukl Il995t jChakrabarti fc Nandil 
120011 iMandal fc Chakrabartil^OloT . As seen from Ta¬ 
ble [2 the intermediate states have comparable values of 
sub-Keplerian (m/,) and Keplerian (rhd) accretion rates, 
which can be due to the c onversio n of one type o f m atte r 
flow into the other type (jMandal fc ChakrabartilfeOlO ). 
Accordingly, a super-Eddington accretion rate may not 
be a pre-requisite for these variability classes to occur. 


The interplay between the different accreting (Keplerian 
and sub-Keplerian) matter with the outflowing matter 
may itself cause such variabilities. It is possible that 
such an interplay between the accreting and outflowing 
streams is currently happening in GRS 1915+105 and 
that in the past this system too might have undergone 
an evolution from the Hard to the Intermedi ate states to 
its pre sent phase like oth er such sources (seelNandi et all 
(2012) for GX 339-4 and lCapitanio et all (I2012D for IGR 
J17091-3624). 

Details of such similarities between the variability 
classes in GRS 1915+105 and those observed in other 
such sources will be explored later and presented else¬ 
where. Until we know a concrete explanation for the 
occurrence of such variability classes, concluding on the 
nature of sources showing such variabilities is difficult. 
However, our current work suggests the presence of a 
high mass stellar black hole in the binary system IGR 
J17091-3624, which accretes at sub-Eddington rates and 
still shows such variability classes. 
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APPENDIX 

APPENDIX I : FITTING WITH ERRORS IN TWO VARIABLES 

This gives a brief outline of the error in variab l e method, and the t echnique that we have used here. For further details, 
please refer to Macdon ald fc Thompsonl (1992); P ress et ahl (|1992L and references therein). One way to accommodate 
erro rs in t h e second varia ble is to alter the fit statistic used for minimization. Various formulations are discussed in 
iMacdonald fc Thompsonl (|1992H . but we use the one (denoted here as O) and given as 

C 2 y (w/fo)) 2 

^ varfyi - f(xi)) ' 

For estimating confidence intervals on the mass, we do Monte Carlo runs of the same fitting routine with the ob¬ 
servations (xi ,yi ) randomly generated as per their error variance fsee lPress et all[1992L for details). We then use a 
Gaussian kernel density estimator to find the probability density function of the obtained set of mass values. The 
errors on the final set of parameters are quoted as ler (68%) deviations from their central value. 

APPENDIX II : PDF CALCULATION FOR THE METHODS MENTIONED IN §3.1 - §3.3 

This section gives an outline of the methods used to construct a PDF from the results of the statistical fitting 
routines. 

• Construct PDF for §3.1 

The PDF is obtained by doing Monte-Carlo simulation runs on the data. We simulate multiple data-sets (~ 10 5 ) 
from the assumed Gaussian error distributions on photon index (r) and QPO frequency (v) (Appendix I). For 
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Figure 5. Chi-square variation plot of parameter (v) used to obtain the PDF. Plot shows different levels of \' 2 , namely Ai at which the 
confidence levels and the confidence intervals (v\,V 2 ) are estimated. The PDF is calculated from these confidence levels. 

each simulated data-set we repeat the least squares fitting routine and note the parameter and statistic value. 
Finally, we obtain a distribution of our parameter of interest (i.e., Mass) by using a Gaussian Kernel Density 
Estimator. 

• Construct PDF for §3.2 and §3.3 

The PDF is obtained from the variation of the fit-statistic (chi-square) about the parameter minima. This gives 
the confidence intervals (rq,tq) on a single parameter from the least-squares fit as seen in equation below 

fV 2 

C v v l = / p(v)dv = P(x 2 > Ai). 

J V\ 


Here, y 2 denotes the chi-squared distribution with 1 degree of freedom, p{v) denotes the PDF of parameter v, 
C^\ denotes the confidence level and Ai is the level (corresponding to V\ and V 2 ) of variation in the y 2 statistic 
(see Figure [5] for details). 

These confidence levels are taken to be the representative areas under the PDF curve in the given interval (iq 
to V 2 ). We find the half area (denoted as A) under the PDF curve from each side of the interval (tq or tq) to 
the minimum value ( v m i n ) by 


AVmin _ Vmin Vl p,v 2 . av 2 _ Vmin p ,„ 2 


V2 — Vi 


V 2 - Vx 


This is applied to Ai which is closest to v m i n . For subsequent A and v. we find the corresponding half areas ( A ), 
by cumulatively adding to the previous estimate of half area (A). Thus, we bin the A values in small steps to 
estimate these half areas for many points. From these half areas (A), the Cumulative Density Function (CDF) 
for each distribution can be found as, 


F(v) 


p{v)dv 


A v rnin- A ^ i ,, ,<v m in 
\ A Vmdn + A% . ,V>V min ' 


Once, we obtain the CDF, the PDF can be calculated by taking a derivative of this, which can be written as : 


p{v) 



APPENDIX III : COMBINING THE ESTIMATES 

The combination of estimates can be done using the likelihood P(Mi\Xi), which denotes the PDF of mass from 
observations A,; of method i (see Appendix II for details on estimation of the PDF). We follow the approach as 
explained below to find the combined / overall limits of the mass from independent PDF only. Thereby, we combine 
the estimate in §3.1 with §3.3 and §3.2 with §3.3 separately. 

• Joint Likelihood estimation using the Naive Bayes theorem 
The joint likelihood of the three methods for BH mass estimation can be obtained using the simple procedure 
known as the ‘naive Bayes’ algorithm where 

P{M) = K\[P{M i \Xi). 
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The Naive Bayes appro ach has been used commonly for multivariate classification of objects like stars in catalogs 
fsee iBroos et al.l l2013f ). We use it here, for parameter confidence level estimation. Here, P{M) denotes the 
combined / final PDF of mass and K is the normalization constant to get unit area PDF. This method, when 
applied to Gaussian distributions gives exactly the weighted average estimate, where the mean of the joint 
PDF will give the weighted average and the standard deviation of the joint PDF will give the error on the 
weighted average. Thus, in this method the width of the final PDF is always less than that of each individual 
PDF. However, this method can only be applied to PDFs which are independent to each other. In classification 
schemes, usage of t he Naive Bayes approach has shown to give acceptable results even if independence is violated, 
as is shown in iDomingos fc Pazzanil (119971) . In our case, if dependence among our mass estimation methods 
exists and is not accounted for, leads to underestimation of the width of the final PDF. To overcome this, we 
do not combine the PDFs from the methods which may have some level of dependence. Another limitation of 
this method is that it gives an offset / biased PDF, if any of the individual PDFs have unaccounted biases / 
systematics. To calculate the joint PDF we perform a multiplication of the individual PDFs for every value of 
mass in the range from 6 M 0 to 20 Mg. 


Once we obtain the final PDF P(M), we quote the confidence interval of the mass by finding the limits beyond 
which the PDF encompasses 5% of the distribution on each side, i.e., we compute Mi and M 2 to give P{M < Mi) = 
P(M > M 2 ) = 0.05. Thus, we find the 90% confidence interval on the mass. We note that this interval may extend 
to a smaller level of confidence than 90% due to an increase in the width of the combined PDF in case of unaccounted 


dependence between the individual PDFs. 
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